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These lecture notes provide an introduction to quantum cluster methods for strongly correlated 
systems. Cluster Perturbation Theory (CPT), the Variational Cluster Approximation (VCA) and 
Cellular Dynamical Mean Field Theory (CDMFT) are described, as well as the exact diagonalization 
solver for the cluster. Potthoff's self-energy functional formalism is reviewed. Some numerical 
procedures, in particular regarding the exact diagonalization method and the frequency-momentum 
integrals needed in VCA, are discussed in detail. 



I. INTRODUCTION 

Classic numerical approaches to lattice models such as 
the Hubbard model are usually based on a solution of the 
model on a periodic lattice with a small number of sites. 
For instance, Exact Diagonalizations (ED) calculations 
are performed on periodic systems with no more than 
<~ 20 sites and Quantum Monte Carlo (QMC) calcula- 
tions are limited in practice on systems with < 100 sites. 
Some extrapolation is then needed to make statements 
about the thermodynamic (e.g. infinite-size) limit. One 
advantage of such approaches is their relative simplicity 
and lack of ambiguity. A disadvantage is that broken 
symmetry states need careful analysis to be identified, as 
they are fully revealed only in the thermodynamic limit. 

Quantum cluster methods are a set of closely related 
approaches that instead consider a finite cluster of sites 
embedded in the infinite lattice. The embedding is done 
by adding to the cluster additional fields or bath degrees 
of freedom such as to best represent the effect of the sur- 
rounding infinite lattice. Variational or self-consistency 
principles are used to set the values of these extra param- 
eters. In these approaches, broken symmetry states can 
appear even for the smallest clusters used, somewhat like 
ordinary mean- field theory. However, unlike mean field 
theory, these approaches are dynamical and retain the 
full effect of strong correlations. These methods are usu- 
ally known by their acronyms: 

1. VCA (Variational Cluster Approximation) 1 or VCPT 
(Variational Cluster Perturbation Theory) 

2. CDMFT (Cluster/Cellular Dynamical Mean-Field 
Theory) 2 

3. DCA (Dynamical Cluster Approximation) 3 ' 4 

The first two of these methods (VCA and CDMFT) can 
be understood within a more general framework called 
the Self-Energy Functional Approach (SFA) , proposed by 
M. Potthoff. 5 The last one (DCA) cannot, as usually for- 
mulated, but is a momentum-space analog of CDMFT. 



Both DCA and CDMFT are cluster generalizations of 
Dynamical Mean Field Theory (DMFT). 6 ' 7 

These lecture notes will concentrate on VCA, its pre- 
cursor CPT, and CDMFT, along with the exact diagonal- 
ization solver. DCA will be discussed by M. Jarrel, later 
during this school. Readers are referred to the excellent 
review by Maier et al. 8 for different aspects of cluster 
methods, including alternate solvers (that review, how- 
ever, was written before the VCA technique was mature). 

Each of these cluster methods is in turn dependent on 
a solution of the cluster Hamiltonian H' which differs 
from the lattice Hamiltonian H by a number of differ- 
ent (exact or approximate) methods. The cluster being 
often compared to an impurity, we often refer to these as 
different impurity solvers, although the expression clus- 
ter solver is more appropriate. In these notes, we will 
describe exact diagonalizations as an cluster solver. 

We will be concerned with the one-band Hubbard 
model, defined on a lattice 7 whose sites will be labelled 
by latin indices {i, j, . . .). The destruction operator for 
an electron on the Wannier orbital centered at the site 
Yi with spin a will be denoted ci a , and the correspond- 
ing number operator will be ri{ a . With this notation, the 
lattice Hamiltonian reads 

H = Y ***4t c j> + u Y n ^ n a - v Y n * w 

i,j,(7 i i 

where tij is the hopping matrix, U is the one-site 
Coulomb repulsion and fi is the chemical potential, which 
we find convenient to include in the Hamiltonian. We will 
assume, for counting purposes, that the lattice 7 is peri- 
odic, with a large (i.e., billions) but finite number of sites 
N. 

This paper is organized as follows: 

1. Section 2 reviews Cluster Perturbation Theory (CPT), 
the simplest of all cluster approaches, which is the 
basis of VCA and serves as a general introduction to 
cluster kinematics. 

2. Section 3 reviews the Exact Diagonalization technique 
for computing the cluster's ground state and Green 
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FIG. 1: A 10-site cluster and the corresponding superlattice 
vectors. 



function, making use of the Lanczos and Band Lanczos 
methods. 

3. Section 4 reviews Potthoff 's self-energy functional ap- 
proach, necessary to understand VCA (Section 5) and 
CDMFT (Section 6). 

4. Roughly a third of this paper consists of appendices 
that explain some specific points in more detail. In 
particular, Appendix A deals with cluster kinematics 
and is required reading before Section 2. 

II. CLUSTER PERTURBATION THEORY 

The simplest quantum cluster method is Cluster 
Perturbation Theory (CPT). 9 ' 10 CPT can be viewed 
as a cluster extension of strong-coupling perturbation 
theory 11 , although limited to lowest order. 12 Its kinemat- 
ical features are found in more sophisticated approaches 
like VCA or CDMFT. The reader is strongly encouraged 
to read Appendix A, where much of the notation about 
clusters and indices is explained. 

CPT proceeds as follows. First a cluster tiling is chosen 
(sec, e.g., Fig. 1). Then the lattice Hamiltonian H is writ- 
ten as H = H' + V, where H' is the cluster Hamiltonian, 
obtained by severing the hopping terms between differ- 
ent clusters, and V contains precisely those terms. V is 
treated as a perturbation. It can be shown, by the tech- 
niques of strong-coupling perturbation theory 10,12 , that 
the-lowest order result for the lattice Green function is 

G^M = G' _1 (w) - V , (2) 

where V is the matrix of inter-cluster hopping terms and 
G'(cj) the exact Green function of the cluster. This for- 
mula deserves a more thorough description: G, G' and 
V are matrices in the space of sites indices i. In terms 
of compound cluster/cluster-site indices (m, a), G' is di- 
agonal in m and identical for all clusters, whereas V is 



essentially off-diagonal in m. Because of translation in- 
variance on the superlattice, the above formula is sim- 
pler in terms of reduced wavevectors, following a partial 
Fourier transform m — * k: 

G- 1 (k,w) = G'- 1 (w)-V(k) . (3) 

The matrices appearing in the above formula are now of 
order L (the number of sites in the cluster), i.e., they 
are matrices in cluster site indices (a, 6, ... ) only (let us 
say we ignore spin or possible band indices for the mo- 
ment). G' is independent of k, whereas V is frequency 
independent. 

The basic CPT relation (3) may also be expressed in 
terms of the self-energy Z of the cluster Hamiltonian as 

G- 1 (k,w) = G - 1 (k,a;)-ZM , (4) 

where Grj(k, oj) is the Green function associated with the 
non-interacting part of the lattice Hamiltonian. This fol- 
lows simply from the relations 

G'- 1 = w - 1' - Z (5) 
Go" 1 = u - t' - V , (6) 

where t' is the restriction to the cluster of the hopping 
matrix (chemical potential included). It is in the form 
(4) that CPT was first proposed 9 . 

A supplemental ingredient to CPT is the periodization 
prescription, that provides a fully k-dependent Green 
function out of the mixed representation G a &(k,aj). In- 
deed, the cluster decomposition breaks the original lat- 
tice translation symmetry of the model. The Green 
function (3) is not fully translation invariant. This 
means that it is not diagonal when expressed in terms 
of wavevectors: G — > G(k,k'). Due to the residual su- 
perlattice translation invariance, however, k' and k must 
differ by an element of the reciprocal superlattice. The 
periodization procedure proposed in Ref. 10 applies to 
the Green function itself: 

G cpt (k, oj) = i Y, e- lk - {r "- rb) G ab (k, u) . (7) 

a.b 

Moreover, since the reduced zone k is taken from is im- 
material, on may replace k by k in the above formula 
(i.e. replacing k by k + K yields the same result). This 
periodization formula may be heuristically justified as fol- 
lows. In the (K,k) basis, the matrix G has the following 
form (see. Eqs (A5)): 

Gkk' (k, lo) = ^e- 4 ( K - r "- K '- r ^G Q ,(k, W ) . (8) 

a, b 

This form can be further converted to the full wavevector 
basis (k = K + k) by use of the unitary matrix A of 
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Eq (A20): 

G(k + K, k + K') = (A(k)GA f (k) 



KK' 



L 2 



g-iti+K-KOT.gitk+K'-K'jTi,^^, 



o,6,Ki,Ki 



= i X) e ~ <c£+K) " ae<c£+K ' ) ' rfcG «*( fc ' w ) 



(9) 



a . b 



The periodization prescription (7) amounts to picking 
the diagonal piece of the Green function (k = k') and 
discarding the rest. This makes sense in as much as the 
density of states N(u>) is the trace of the imaginary part 
of the Green function: 

N(w) = -21m tr G(w) = -21m ^ G lt {uj) 

i 

= -2Im^G(k, W ) , (10) 

k 

and the spectral function A(k,uj), as a partial trace, in- 
volves only the diagonal part. Indeed, it is a simple mat- 
ter to show from the anticommutation relations that the 
frequency integral of the Green function is the unit ma- 
trix: 



-21m 



2^ 



G(w) = 1 



(11) 



This being representation independent, it follows that 
the frequency integral of the imaginary part of the off- 
diagonal components of the Green function vanishes. 

Another possible prescription for periodization is to 
apply the above procedure to the self-energy X instead. 
This is appealing since H is an irreducible quantity, as 
opposed to G. This amounts to throwing out the off- 
diagonal components of X before applying Dyson's equa- 
tion to get G, as opposed to discarding the off-diagonal 
part at the last step, once the matrix inversion towards G 
has taken place. As Fig. 2 shows, periodizing the Green 
function (Eq. (7)) reproduces the expected feature of the 
spectral function of the one-dimensional Hubbard model. 
In particular, the Mott gap that opens at arbitrary small 
U (as known from the exact solution), whereas periodiz- 
ing the self-energy leaves spectral weight within the Mott 
gap for abritrary large value of U. Therefore we will al- 
ways use Green function periodization. 

CPT has the following characteristics: 

1. Although it is derived using strong-coupling pertur- 
bation theory, it is exact in the £/ — > limit, as the 
self-energy disappear in that case. 

2. It is also exact in the strong-coupling limit Uj/U — > 0. 

3. It provides an approximate lattice Green function for 
arbitrary wave vectors. Hence its usefulness in com- 
paring with ARPES data. Even though CPT does 
not have the self-consistency present in DMFT type 




FIG. 2: Top: CPT spectral function of the one-dimensional, 
half-filled Hubbard model with U = 4, t = 1, with Green 
function periodization (L — 16). Bottom : the same, with 
Self-energy periodization instead; notice the important spec- 
tral weight in the middle of the Mott gap. 



approaches, at fixed computing resources it allows 
for the best momentum resolution. This is particu- 
larly important for the ARPES pseudogap in electron- 
doped cuprates that has quite a detailed momentum 
space structure, and for d-wave superconducting cor- 
relations where the zero temperature pair correlation 
length may extend well beyond near-neighbor sites. 

Although formulated as a lowest-order result of strong- 
coupling perturbation theory, it is not controlled by 
including higher-order terms in that perturbation ex- 
pansion - this would be extremely difficult - but rather 
by increasing the cluster size. 

It cannot describe broken-symmetry states. This is 
accomplished by VCA and CDMFT, which can both 
be viewed as extensions or refinements of CPT. 



III. EXACT DIAGONALIZATIONS 

Before going on to describe more sophisticated quan- 
tum cluster approaches, let us describe in some detail a 
particular cluster solver, i.e., a particular method used 
to calculate the ground state and Green function of the 
cluster: the exact diagonalization method, based on the 
Lanczos algorithm. The quantum cluster methods de- 
scribed here are not tied to a specific solver for the clus- 
ter. For instance, Quantum Monte Carlo or any other 
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approximate method of solution for the cluster Green 
function could be used. The exact diagonalization (ED) 
method has the advantage of high numerical accuracy at 
zero temperature, and can be to some extent controlled 
by the size of the cluster used. 

The basic idea behind exact diagonalization is one of 
brute force, but its practical implementation may require 
a lot of care depending on the desired level of optimiza- 
tion. Basically, an exact representation of the Hamil- 
tonian action on arbitrary state vectors must be coded 
- this may or may not involve an explicit construction 
of the Hamiltonian matrix. Then the ground state is 
found in an quasi-exact way by an iterative method such 
as the Lanczos algorithm. The Green function is there- 
after calculated by similar means to be described below. 
The main difficulty with execution is the large memory 
needed by the method, which grows exponentially with 
the number of degrees of freedom. As for coding, the 
main difficulty is to optimize the method, in particular 
by taking point group symmetries into account. 



A. Coding of the basis states 

The first step in the exact diagonalization procedure is 
to define a coding scheme for the quantum basis states. 
A basis state may be specified by the occupation number 
n p (= or 1) of electrons in the orbital labelled \i and has 
the following expression in terms of creation operators: 

(ci T )^ . . . (c^r^t^u . . . (c t. i)niJ | 0) (12) 

where the order in which the creation operators are ap- 
plied is a matter of convention, but important. If the 
number of orbitals is smaller than or equal to 32, the 
string of n M 's forms the binary representation of a 32-bit 
unsigned integer b, which can be split into spin up and 
spin down parts: 



b = b ] + 2 L b i 



(13) 



There are 2 2L such states, but not all are relevant, since 
the Hubbard Hamiltonian is block-diagonal : The num- 
ber of electrons of a given spin (iVj and N± ) is conserved 
and commutes with the Hamiltonian H . Therefore the 
exact diagonalization is to be performed in a sector (i.e. 
a subspace) of the total Hilbert space with fixed values of 
iVf and N± . This space has the tensor product structure 



V = V N , ® v Nl 



and has dimension d — d(N^)d(Ni), where 

L\ 

d{N a ) 



N a \{L-N a )\ 



(14) 



(15) 



is the dimension of each factor, i.e., the number of ways 
to distribute N a electrons among L sites. 



Note that the ground state |f2) of the Hamiltonian 
generally belongs to the sector = N^. For a half- 
filled, zero spin system (iVj = N± = L/2), this trans- 
lates into d = (L!/(i/2)! 2 ) 2 , which behaves like 4 L /L for 
large L: The size of the eigcnproblcm grows exponentially 
with system size. By contrast, the non-interacting prob- 
lem can be solved only by concentrating on one-electron 
states. For this reason, exact diagonalization of the Hub- 
bard Hamiltonian is restricted to systems of the order of 
16 sites or less. 

In practice, a generic state vector is represented by an 
(i-component array of double precision numbers. In or- 
der to apply or construct the Hamiltonian acting on such 
vectors, we need a way to translate the label of a basis 
state (an integer i from to d— 1), into the binary repre- 
sentation (12). The way to do this depends on the level 
of complexity of the Hilbert space structure. In the sim- 
ple case (14), one needs, for each spin, to build a two-way 
look-up table that tabulates the correspondence between 
consecutive integer labels and the binary representation 
of the spin up (resp. spin down) part of the basis state. 
Thus, given a binary representation (6-p, 6^) of a basis 
state | b) = one immediately finds integer labels 

7j(6j) and and the label of the full basis state may 

be taken as 



7 T (6 T ) +d N ^I l (b l ) 



(16) 



On the other hand, given a label i, the corresponding 
labels of each spin part are 



mod (i, djvf) ii—i/d N i : 



(17) 



where integer division (i.e. without fractional remainder) 
is used in the above expression. The binary representa- 
tion b is recovered by inverse tables B as 



B T (i T ) 



(18) 



The next step is to construct the Hamiltonian matrix. 
The particular structure of the Hubbard model Hamil- 
tonian brings a considerable simplification in the simple 
case studied here. Indeed, the Hamiltonian has the form 



H = X T ® 1 + 1 ® K i + V in 



(19) 



where only acts on up electrons and on down elec- 
trons, and where the Coulomb repulsion term Vj nt . is di- 
agonal in the occupation number basis. Thus, storing the 
Hamiltonian in memory is not a problem : the diagonal 
Vi n t. is stored (an array of size d), and the kinetic energy 
K a (a matrix having a small fraction of d\ elements) is 
stored in sparse form. Constructing this matrix, formally 
expressed as 



K 



*ai>4 C & 

a.b 



(20) 



needs some care with the signs. Basically, two basis states 
\b a ) and \b' a ) are connected with this matrix if their bi- 
nary representations differ at two positions a and b. The 
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matrix clement is then (— l) Mab t a i,, where M a f, is the 
number of occupied sites between a and b, i.e., assum- 
ing a < b, 



6-1 



M ab = J2 



(21) 



c=a-\-l 



For instance, the two states (10010110) and (10011100) 
with L = 8 are connected with the matrix element +t^, 
where the sites are numbered from to L — 1. 

Calculating the Hubbard interaction is straighforward: 
a bit-wise AND is applied to the up and down parts of a 
binary state (6| & b^ in C or C++) and the number of set 
bits of the result is the number of doubly occupied sites 
in that basis state. 



B. The Lanczos algorithm for the ground state 

Next, one must apply the exact diagonalization 
method per se, using the Lanczos algorithm. Generally, 
the Lanczos method 13 is used when one needs the ex- 
treme eigenvalues of a matrix too large to be fully di- 
agonalized (e.g. with the Householder algorithm). The 
method is iterative and involves only multiply-add's from 
the matrix. This means in particular that the matrix 
does not necessarily have to be constructed explicitly, 
since only its action on a vector is needed. Thus, is some 
extreme cases where it is practical to do so, the matrix 
elements can be calculated 'on the fly', and this allows 
to save the memory associated with storing the matrix 
itself. 

The basic idea behing the Lanczos method is to build 
a projection J^f of the full Hamiltonian matrix H onto 
the so-called Krylov subspace. Starting with a (random) 
state \4>o), the Krylov subspace is spanned by the iterated 
application of H : 

,X-span{|0 o ),i/|0 o ),i/ 2 |0 o ),... ,H M \4> Q )} (22) 

the generating vectors above are not mutually orthogo- 
nal, but a sequence of mutually orthogonal vectors can 
be built from the following recursion relation 

|0 n+ i) = H\4> n ) - a n \<t> n ) - l£|^ n _i) (23) 



where 



{<t>n\<j>n) 

{<t>n\(i>n) "™ (</>„_! |</>„_l) 



b n \H\(j) n ) b2 



bo = (24) 



and we set the initial conditions bo = 0, \<f>—i) = 0. At 
any given step, only three state vectors are kept in mem- 
ory (</>„ + i, (f> n and (f) n -i). In the basis of normalized 
states \n) — \<j> n )/y/ {4>n\4>n), the projected Hamiltonian 



has the tridiagonal form 

(a b 1 

h ai b 2 

H = b 2 a 2 b 3 



\0 








(25) 



a N J 



Such a matrix is readily diagonalized by fast methods 
dedicated to tridiagonal matrices, and a convergence cri- 
terion must be set for the lowest eigenvalue Eq, at which 
iterations stop. For instance, one may stop the procedure 
when the lowest eigenvalue Eq changes by no more that 
one part in 10 12 . This may require between a number M 
of iterations between a few tens and ~ 200, depending 
on system size. 

The ground state energy E and the ground state \Q) 
are very well approximated by the lowest eigenvalue and 
the corresponding eigenvector of that matrix, which are 
obtained by standard methods. This provides us with 
the ground state |f2) in the reduced basis {| </>«)}• But 
we need the ground state in the original basis, and this 
requires retracing the Lanczos iterations a second time - 
for the \4> n ) are not stored in memory - and constructing 
the ground state progressively at each iteration from the 
known coefficients (f2| </>„). 

The Lanczos procedure is simple and efficient. Conver- 
gence is fast if the lowest eigenvalue E is well separated 
from the next one (Ei). It slows down if E\ — Eq is 
small. If the ground state is degenerate (E\ = E ), the 
procedure will converge to a vector of the ground state 
subspace, a different one each time the initial state \4>q) 
is changed. 

Note that the sequence of Lanczos vectors \<f> n ) is in 
principle orthogonal, as this is garanteed by the three- 
way recursion relation (23). However, numerical error 
will introduce 'orthogonality leaks', and after a few tens 
of iterations the Lanczos basis will become overcomplete 
in the Krylov subspace. This will translate in multiple 
copies of the ground state eigenvalue in the tridiagonal 
matrix (25), which should not be taken as a true degen- 
eracy. However, as long as one is only interested in the 
ground state and not in the multiplicity of the lowest 
eigenvalues, this is not a problem. 



C. The Lanczos algorithm for the Green function 

Once the ground state is known, it remains to calcu- 
late the cluster Green function. The zero-temperature 
Green function G^ v (u}) has the following expression, as 
a function of the complex- valued frequency u>: 



G^M = (Q\d 



UJ - H+ E 

1 

W + H -E 



ct\n) 

c„|ft> 



(26) 
(27) 

(28) 



G 



In the basic Hubbard model, spin is conserved and we 
need only to consider the creation and annihilation of 
up-spin electrons. 

We will first describe a Lanczos algorithm for calcu- 
lating the Green function, that provides a continued- 
fraction representation of its frequency dependence. In 
the next subsection, we will instead present an alternate 
method based on the Band Lanczos algorithm, that pro- 
vides a Lehmann representation of the Green function 
and that is both faster and more memory intensive. 

Consider first the function G^ e (w). One needs to 
know the action of (u> — H + i?o) _1 on the state = 
cjjfi), and then to calculate 

As with any generic function of H , this one can be ex- 
panded in powers of H: 

-^- = - + \h+^H 2 + --- (30) 

z — H z z z z i 

and the action of this operator can be evaluated ex- 
actly at order H M in a Krylov subspace (22). Thus 
we again resort to the Lanczos algorithm: A Lanczos 
sequence is calculated from the initial, normalized state 
\4>a) = |<?V)/&o where bl = (^|^). This sequence gen- 
erates a tridiagonal representation of H , albeit in a dif- 
ferent Hilbert space sector : that with jVj + 1 up-spin 
electrons and N± down-spin electrons. Once the preset 
maximum number of Lanczos steps, or a near zero value 
of b n , has been reached, the tridiagonal representation 
(25) may then be used to calculate (29). This amounts 
to the matrix element Oq[(uj ~ H + £7o) _1 ]oo (the first el- 
ement of the inverse of a tridiagonal matrix), which has 
a simple continued fraction form : 14 

G'^M = ^2 (31) 

0,-ao 

lu — a\ 

U — 0,2 

Thus, evaluating the Green function, once the arrays 
{a n } and {b n } have been found, reduces to the calcu- 
lation of a truncated continued fraction, which can be 
done recursively in M steps, starting from the bottom 
floor of the fraction. 

Consider next the case fi^f. The continued fraction 
representation applies only to the case where the same 
state \4>) appears on the two sides of (29). If fx ^ v, this 
is no longer the case, but we may use the following trick : 
we define the combination 

G+ » - (n\(c, + c,) uj _^ + Eq {c, + c^\Q) (32) 
Using the symmetry G^ Vie (ui) = G„ M , e (w), this leads to 
GWM - ^(G+ » - G^H - <?„„,») (33) 



where G+, e can be calculated in the same way as G MAl>e , 
i.e., with a simple continued fraction. We proceed like- 
wise for G+„ jfc (w). 

Thus, the cluster Green function is encoded in L(L + 
1) continued fractions, whose coefficients are stored in 
memory, so that G'(w) can be computed on demand for 
any complex frequency u>. 

Note that a minimal way to take advantage of clus- 
ter symmetries is to restrict the calculation of the Green 
function to an irreducible set of pairs (fi, v) of orbitals 
that can generate all other pairs by symmetry operations 
of the cluster. Thus, if a symmetry operation g takes the 
orbital \x into the orbital g(jJL), we have 

G'^)=G' gWg{v) {uj) (34) 

Taking this into account is an easy and important time 
saver. 



D. The Band Lanczos algorithm for the Green 
function 

An alternate way of calculating the cluster Green func- 
tion is to apply the band Lanczos procedure 15 . This is 
a generalization of the Lanczos procedure in which the 
Krylov subspace is spanned not by one, but by many 
states. Let us assume that up and down spins are decou- 
pled, so that the Green function is L x L block diagonal. 
The L states = c* |fi) are first constructed, and then 
one builds the projection of H' on the Krylov sub- 
space spanned by 

fl^),...,!^),^'!^),...,^'!^),..., 

1 , (35) 

{H>) M \<t>i),-..,{H') M \H)} 

A Lanczos basis {\n}} is constructed by successive appli- 
cation of H' and orthonormalization with respect to the 
previous 2L basis vectors. In principle, each new basis 
vector \n) is already automatically orthogonal to basis 
vectors |1) through \n— 2L— 1), although 'orthogonality 
leaks' arise eventually and may be problematic. A prac- 
tical rule of thumb to avoid these problems is to control 
the number M of iterations by the convergence of the 
lowest eigenvalue of (e.g. to one part in 10 10 ). In- 
dependently of this, one must be careful about potential 
redundant basis vectors in the Krylov subspace, which 
must be properly 'deflated'. 15 The number of states R in 
the Krylov subspace at convergence is typically between 
100 and 300, depending on system size. The Rx R ma- 
trix J$?, which has a tridiagonal structure in the ordinary 
Lanczos method, now has a band structure made of 2L 
diagonals around the central diagonal. It is then a simple 
matter to obtain a Lehmann representation of the Green 
function in the Krylov subspace (see Appendix B) by cal- 
culating the projections Q^ r of on the eigenstates 
of (the inner products of the |0u)'s with the Lanc- 
zos vectors are calculated as the latter are constructed). 
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The Green function can then be expressed in a Lehmann 
representation (B4) . The two contributions G'^ v e and 
G' h to the Green function are computed separately, 
and the corresponding matrices Q and A are simply con- 
catenated to form the complete Q- and A-matrices, which 
are then stored and allow again for a quick calculation 
of the Green function as a function of the complex fre- 
quency w. The matrix 2L x R matrix Q has the property 
that 

QQ f = 12LX2L (36) 

This holds even if the Lehmann representation is ob- 
tained from a subspace and not the full space, and is 
simply a consequence of the anticommutation relations 

The band Lanczos method requires more memory than 
the usual Lanczos method, since 2L + 1 vectors must si- 
multaneously be kept in memory, compared to 3 for the 
simple Lanczos method. On the other hand, it is faster 
since all pairs (n, v) are covered in a single procedure, 
compared to L(L + l)/2. Thus, we gain a factor I? in 
speed at the cost of a factor L in memory. Another ad- 
vantage is that it provides a Lehmann representation of 
the Green function. 



E. Cluster symmetries 

It is possible to optimize the exact diagonalization pro- 
cedure by taking advantage of the symmetries of the 
cluster Hamiltonian, in particular coming from cluster 
geometry. If the Hamiltonian is invariant under a dis- 
crete group <S of symmetry operations and |(S| denotes 
the number of such elements (the order of the group), 
the dimension of the largest Hilbert space needed can 
be reduced by a factor of almost |25|, and the number 
of state vectors needed in the band Lanczos method re- 
duced by the same factor. The corresponding speed gain 
is appreciable. In the case of large clusters (e.g. 16 sites), 
taking advantage of symmetries may make the difference 
between doing or not doing the problem. The price to 
pay is a higher complexity in coding the basis states, 
which almost forces one to store the Hamiltonian ma- 
trix in memory, if it were not already, since calculating 
matrix elements 'on the fly' becomes more time consum- 
ing. Note that we are using open boundary conditions 
(except in the case of the DCA, not discussed in these 
notes), and therefore there is no translation symmetry 
within the cluster; thus we are concerned with points 
groups, not space groups. 

Let us start with a simple example: a cluster invari- 
ant with respect to a single inversion, or a single rotation 
by 7r. One may think of a one-dimensional cluster, for 
instance, with a left-right inversion. The corresponding 
symmetry group is C2, with two elements: the identity 
e and the inversion 1. The group C2 contains two irre- 
ducible representations, noted A and B, corresponding 
respectively to states that are even and odd with respect 



to l. Because the Hamiltonian is invariant under inver- 
sion: H = l Hl, eigenvectors of H will be either even 
or odd, i.e. belong either to the A or to the B represen- 
tation. Likewise, the Hamiltonian will have no matrix 
elements between states belonging to different represen- 
tations (the reader is invited to read Appendix C for a 
review of the necessary group-theoretical concepts) . 

In order to take advantage of this fact, one needs to 
construct a basis containing only states of a given repre- 
sentation. The occupation number basis states |6) (or bi- 
nary states, as we will call them) introduced above are no 
longer adequate. In the case of the simple group C2, one 
should rather consider the even and odd combinations 
I b) ± b\b) (and some of these combinations may vanish). 
Yet we still need a scheme to label the different basis 
states and have a quick access to their occupation num- 
ber representation, which allows us to compute matrix 
elements. Let us briefly describe how this can be done (a 
more detailed discussion can be found, e.g., in Ref. 16). 
Under the action of the group ©, each binary state gener- 
ates an 'orbit' of binary states, whose length is the order 
|<S| of the group, or a divisor thereof. To such an orbit 
corresponds at most d a states in the irreducible represen- 
tation labeled a, given by the corresponding projection 
operator: 

i^ = ||E4 a) *#> (37) 
1 1 9 

where d a is the dimension of the irreducible representa- 
tion a. We will restrict the discussion to the simplest 
case, where all irreducible representations considered are 
one-dimensional (d a — 1; the case d a > 1 turns out to be 
quite a bit more complex). Then the state \\p) is either 
zero or unique for a given orbit. We can then select a 
representative binary state for each orbit (e.g. the one 
associated with the smallest binary representation) and 
use it as a label for the state \ip). We still need an index 
function B(i) which provides the representative binary 
state for each consecutive label i. The reverse correspon- 
dence i = 1(b) is trickier, since symmetrized states are no 
longer factorized as products of up and down spin parts. 
It is better then to search the array B for the value of the 
index i that provides a given binary state b. One can still 
be aided by a partial reverse index If (6f) that provides 
the first occurence in the list B of a state with b^ as the 
spin up part, assuming that states are sorted according 
to 6|, then according to b±. 

Once the basis has been constructed, one needs to con- 
struct a matrix representation of the Hamiltonian in that 
representation. Given two states \ipi) and \4>2), repre- 
sented by the binary states |6i) and I&2), it is a simple 
matter to show that the matrix element is 

<V>2|i^i) = ^^X^MVWHlh) (38) 

where the phase <fr g (b) is defined by the relation 

g\b) = Mb)\9b) ■ (39) 
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TABLE I: Number of matrix elements of a given value in the 
nearest-neighbor hopping operator on the half-filled 3x4 = 12 
site cluster, for each irreducible representation of Ci v . The 
dimension of each subspace is indicated on the second row. 







A 2 


Si 


B 2 


dim. 


213,840 


213, 248 


213,440 


213,248 


value 










-2 


96 


736 


704 





-V2 


12,640 


6,208 


7, 584 


5,072 


-1 


2, 983, 264 


2,936,144 


2, 884, 832 


2,911,920 


1 


952,000 


997, 168 


1,050,432 


1,021,392 


V2 


5,088 


2,304 


3,232 


2,992 


2 


32 












In the above relation, \gb) is the binary state obtained 
by applying the symmetry operation g to the occupation 
numbers forming b, whereas the phase <fi g (b) is the prod- 
uct of signs collected from all the permutations of cre- 
ation operators needed to go from b to gb. Formula (38) 
is used as follows to construct the Hamiltonian matrix: 
First, the Hamiltonian can be written as H = ^ r H r , 
where H r is a hopping term between specific sites, or a 
diagonal term like the interaction. One then loops over 
all bis. For each b\, and each term H r , one construct 
the single binary state H r \b\). One then finds the rep- 
resentative b 2 of that binary state, by applying on it all 
possible symmetry operations until g is found such that 
\gb2) = H r \b\). During this operation, the phase (f> g (b) 
must also be collected. Then the matrix element (38) is 
added to the list of stored matrix elements. Since each 
term H r individually is not invariant under the group, 
there will be more matrix elements generated than there 
should be, i.e., there will be cancellations between differ- 
ent matrix elements associated with the same pair (pi, 
b 2 ) and produced by the different H r 's. For this reason, 
it is useful to first store all matrix elements associated 
with a given bi in an intermediate location in order for 
the cancellations to take effect, and then to store the 
cleaned up 'column' labelled by b\ to its definitive stor- 
age location. Needless to say, one should only store the 
row and column indices of each clement of a given value. 

Table I gives the values and number of matrix elements 
found for the nearest-neighbor hopping terms on the half- 
filled 12-site (3 x 4) cluster, in each of the four irreducibe 
representations of the group C 2v ■ 

F. Green functions using cluster symmetries 

Most of the time, the ground state lies in the trivial 
(symmetric) representation. However, taking advantage 
of symmetries in the calculation of the Green function re- 
quires all the irreducible representations to be included 
in the calculation. Consider for instance the simple ex- 
ample of a C 2 symmetry, with a ground state |f2) in the 
A (even) representation. Constructing the Green func- 



tion involves applying on |f2) the destruction operator c a 
(or the creation operator c\) associated to site a. The 
excited state thus produced does not belong to a well- 
defined representation. Instead, on should destroy (or 
create) and electron in an odd or even state, by using the 
linear combinations c a ±c ta , where ta is the site obtained 
by applying the symmetry operation to a. Thus, in calcu- 
lating the Green function (26), one should express each 
creation/destruction operator in terms of symmetrized 
combinations, e.g., 

Ca = ^(c a + C M ) + ^(c a - C La ) (40) 

More generally, one would use symmetrized combinations 
of operators 

4 a) = E M P« )c « ( 41 ) 

a 

such that Cp Q ^ transforms under representation a, and p 
labels the different possibilities. For instance, for a linear 
cluster of length 4 and an inversion symmetry that maps 
the sites (1234) into (4321), these operators are 

(A) , (B) 

C] — Cl + Ca C] ' — Cl — Ca 
1 14 1 14 

(A) . (B) 1 1 

c 2 = c 2 + c 3 c\ ' = c 2 ~c 3 

Then, for each representation, one may use the Band 
Lanczos procedure and obtain a Lehmann representa- 
tion Qp" for the associated Green function G P a(u>). If 
the ground state is in representation a and the operators 
c p ^ of representation j3 are used, the Hilbert space sec- 
tor to work with will be the tensor product representation 
a ® (3, which poses no problem at all when all irreps are 
one-dimensional, but would bring additional complexity 
if the ground state were in a multi-dimensional repre- 
sentation. Finally, one may bring together the different 
pieces, by building alxL matrix M pa that is the vertical 
concatenation of the various rectangular matrices Mp%} , 
and returning to the usual Q-matrix representation 

Qar = (M-^apQpr (43) 

Using cluster symmetries for the Green function saves 
a factor |<S| in memory because of the reduction of the 
Hilbert space dimension, and an additional factor of \<3\ 
since the number of input vectors in the band Lanczos 
procedure is also divided by \(&\. Typically then, most of 
the memory will be used to store the Hamiltonian matrix. 

G. Parallelization 

For larger clusters (e.g. 16 sites), the computer mem- 
ory required to carry out the exact diagonalization is too 
large to fit on a standard computer. In those cases the 
only practical choice is to parallelize the program. Al- 
though this is a technical issue that has more to do with 
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FIG. 3: How to split a matrix- vector multiplication \y) = 
H\x) across two processes. 'Blue' data reside on one process, 
and 'red' on the other. For each of the two segments of \y), 
each process performs a block matrix multiplication, and the 
results of the two processes must be transfered to each other 
to be added. 



programming than with the algorithm, a brief explana- 
tion is in order. Parallelization consists in dividing the 
task and data between many processes (run on differ- 
ent cpus) , with communication between processes taking 
place on a frequent basis. The Message Passing Interface 
(MPI) Library is the most common way to accomplish 
this on distributed-memory machines. Parallelization is 
often a difficult task, and is likely not to scale well (i.e., 
the inverse computing time grows more slowly than the 
number of processes) when inter-process communications 
occur too frequently. However, parallelization makes the 
difference between doing or not doing a large problem. 

Let us now briefly describe a possible way to par- 
allelized an exact diagonalization program, as used by 
us. Let ,jV be the number of processes across which the 
problem is parallelized. We split each Hilbert space vec- 
tor into jV (nearly) equal segments, and the Hamilto- 
nian matrix into jV 2 blocks (labelled Hjj, with I, J — 
A single matrix- vector multiplication \y) — 
H\x) then proceeds, for each process, by jV successive 
operations \y)i = Hjj\x)j, J labelling the different pro- 
cesses and / the successive operations. After each oper- 
ation the resulting vectors \y)i must be sent to process 
/ to be summed in a single segment. This is illustrated 
on Fig. 3 for jV = 2. Thus, each multiply-add oper- 
ator involves ,jV 'broadcast' or 'reduce' operations, in 
MPI jargon. The construction of the Hamiltonian is also 
parallelized, as each process takes care of its own group 
of columns. This constitutes what is called fine-grained 
parallelization: communications are very frequent (many 
calls per matrix- vector multiply add). Consequently, 
scaling is poor and in practice the number of processes 
should be kept to a minimum, just enough to fit the pro- 
gram in memory. 

As a whole, computational scientists will feel an ever 
increasing pressure to use parallel computing, as this will 
become the only way not only to do larger problem, but 
to substantially speed up all problems, because of the 
slowing down of Moore's law and of the ubiquity of cpus 
with an increasing number of cores. 



IV. THE SELF-ENERGY FUNCTIONAL 
APPROACH 

That CPT is incapable of describing broken symme- 
tries is its major drawback. Treating spontaneously bro- 
ken symmetries requires some sort of self-consistent pro- 
cedure, or a variational principle. Ordinary mean-field 
theory does precisely that, but is limited by its discard- 
ing of fluctuations and its uncontrolled character. 

A heuristic way of treating broken symmetry states 
within CPT would be to add to the cluster Hamiltonian 
H' a Weiss field that pushes the system towards some 
predetermined form of order. For instance, the follow- 
ing term, added to the Hamiltonian, would induce Neel 
antifcrromagnetism: 

H' M = MO M = M elQ ' ra Kt - n ai ) (44) 

a 

where Q = (ir, tt) is the antiferromagnetic wavevector. 
What is needed is a procedure to set the value of the 
Weiss parameter M. Adopting a mean- field- like pro- 
cedure (i.e. factorizing the interaction in the correct 
channel and applying a self-consistency condition) would 
bring us exactly back to ordinary mean-field theory: the 
interaction having disappeared, the cluster decomposi- 
tion would be suddenly useless and CPT would provide 
the same result regardless of cluster size. 

The solution to that conundrum is most elegantly pro- 
vided by the self-energy functional approach (SFA), pro- 
posed by Potthoff. 1 This approach also has the merit of 
presenting various cluster schemes from a unified point 
of view. It can also be seen as a special case of the more 
general inversion method 17 , recently reviewed in Ref. 18 
in the context of Density Functional Theory and DMFT. 

To start with, let us introduce a functional f2 t [G] of the 
Green function: 

fit[G] = *[G] - Tr ((G^ - G^G) + Tr ln(-G). (45) 

This means that, given any Green function Gij(uj) one 
can cook up - yet with the usual analytic properties of 
Green functions as a function of frequency - this expres- 
sion yields a number. In the above expression, products 
and powers of Green functions - e.g. in series expan- 
sions like that of the logarithm - are to be understood 
in a functional matrix sense. This means that position i 
and time r, or equivalently, position and frequency, are 
merged into a single index. Accordingly, the symbol Tr 
denotes a functional trace, i.e., it involves not only a sum 
over sites indices, but also over frequencies. The latter 
can be taken as a sum over Matsubara frequencies at fi- 
nite temperature, or as an integral over the imaginary 
frequency axis at zero temperature. 

The Luttinger Ward functional $[G] entering this ex- 
pression is usually denned as the sum of two-particle ir- 
reducible (2PI) diagrams : diagrams that cannot by split 
into disjoint parts by cutting two fermion lines (Fig. 4). 
These are sometimes called skeleton diagrams, although 
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FIG. 4: Diagrammatic definition of the Luttinger-Ward func- 
tional, as a sum over skeleton graphs. 



'two-particle irreducible' is more accurate. A diagram- 
free definition of <1>[G] is also given in Ref. 19. For our 
purposes, what is important is that (1) The functional 
derivative of <&[G] is the self-energy 



<5$[G] 
<5G 



= E 



(46) 



(as defined diagramatically) and (2) it is a universal func- 
tional of G in the following sense: whatever the form of 
the one-body Hamiltonian, it depends only on the in- 
teraction and, functionnally, it has the same dependence 
on G. This is manifest from its diagrammatic definition, 
since only the interaction (dotted lines) and the Green 
function given as argument, enter the expression. The 
dependence of the functional fl t [G] on the one-body part 
of the Hamiltonian is denoted by the subscript t and it 
comes only through G^ 1 = w — t appearing on the right- 
hand side of Eq. (45). 

The functional Q t [G] has the important property that 
it is stationary when G takes the value prescribed by 
Dyson's equation. Indeed, given the last two equations, 
the Eulcr equation takes the form 



SQ t [G] 
SG 



= Z - G^ + G- 1 = 0. 



(47) 



This is a dynamic variational principle since it involves 
the frequency appearing in the Green function, in other 
words excited states are involved in the variation. At this 
stationary point, and only there, Sl t [G] is equal to the 
physical (thermodynamic) grand potential. Contrary to 
Ritz's variational principle, this last equation does not 
tell us whether f2 t [G] is a minimum, a maximum, or a 
saddle point there. 

There are various ways to use the stationarity prop- 
erty that we described above. The most common one 
is to approximate $[G] by a finite set of diagrams. 
This is how one obtains the Hartree-Fock, the FLEX 
approximation 20 or other so-called thcrmodynamically 
consistent theories. This is what Potthoff calls a type II 
approximation strategy. 21 A type I approximation sim- 
plifies the Euler equation itself. In a type III approxi- 
mation, one uses the exact form of <1>[G] but only on a 
limited domain of trial Green functions. 

Following Potthoff, we adopt the type III approxima- 
tion on a functional of the self-energy instead of on a 
functional of the Green function. Suppose we can locally 



invert Eq. (46) for the self-energy to write G as a func- 
tional of Z. We can use this result to write, 

n t [Z] = F[Z] - Tr lnt-Go" 1 + Z). (48) 

where we defined 

F[Z] = <j>[G] — Tr(EG). (49) 

and where it is implicit that G = G[E] is now a func- 
tional of Z. F[Z], along with the expression (46) for the 
derivative of the Luttinger-Ward functional, defines the 
Legendre transform of the Luttinger-Ward functional. It 
is easy to verify that 



5F[Z] _ 6$[G] SG[Z] ^SG[Z] 



sz 



SG 5Z 



SZ 



(50) 



hence, flt[Z] is stationary with respect to Z when Dyson's 
equation is satisfied 



SQt[Z] 
SZ 



-G + (G 0t 1 -Z)- 



0. 



(51) 



To perform a type III approximation on F[Z], we take 
advantage that it is universal, i.e., that it depends only 
on the interaction part of the Hamiltonian and not on 
the one-body part. We then consider another Hamilto- 
nian, denoted H' and called the reference system, that 
describes the same degrees of freedom as H and shares 
the same interaction (i.e. two-body) part. Thus H and 
H' differ only by one-body terms. We have in mind for 
H' the cluster Hamiltonian, or rather the sum of all (mu- 
tually decoupled) cluster Hamiltonians. At the physical 
self-energy Z of the cluster, Eq. (48) allows us to write 



n t ,[Z] = Q' = F[Z] - Tr ln(-G') 



(52) 



where il' is the cluster Hamiltonian 's grand potential and 
G' its physical Green function, obtained through the ex- 
act solution. From this we can extract F[Z] and it follows 
that 



n t [Z] = Q'+ Tr ln(-G') 
= SY + Tr ln(-G') 



Tr lnt-Go" 1 + E) 
Tr ln(-G) 



(53) 



where G now stands for the CPT Green function (2). 
This expression can be further simplified as 



fi t [E] = SI' — Tr ln(l - VG') 



(54) 



Let us finally make the trace more explicit: It is a sum 
over frequencies and a sum over lattice sites (and spin and 
band indices), which can be expressed instead as a sum 
over reduced wavevectors (as the CPT Green function is 
diagonal in that index), plus a "small" trace (denoted 
tr) on residual indices (cluster site, spin, and band): 



$l t [E] = Q' tr ln [l- V(k)G'(k,w) 

= S1'-T^^lndct [l-V(k)G'(k,w)l (55) 



w k 



11 



where the matrix identity tr In A = In det A was used in 
the second equation. 

The type III approximation comes from the fact that 
the self-energy Z is restricted to the exact self-energy of 
the cluster problem H' , so that variational parameters 
appear in the definition of the one-body part of H'. To 
come back to the question of the Weiss field M introduced 
at the beginning of this section, we would set its value 
by solving the cluster Hamiltonian - i.e., calculating f2' 
and G' - for many different values of M and evaluate the 
functional (55) for each of them, selecting the value that 
makes Expression (55) stationary. This is the idea behind 
the variational cluster approximation (VCA), described 
in more detail in the next section. 

In practice, we look for values of the cluster one-body 
parameters t' such that <5Q t [5Z]/<5t' = 0. It is useful for 
what follows to write the latter equation formally, al- 
though we do not use it in actual calculations. Given 
that ft' is the actual grand potential evaluated for the 
cluster, dCl' /dt' is canceled by the explicit t' dependence 
of Tr ln(— Gg t , + Z) and we are left with 



= 



Sn t [T.} <5Z 



<SZ 
Tr 



St' 



J 0t' 



Got 1 



(56) 



This may be explicited as 




Got? - E H 



A^VGot 1 ^)-^)^ 



St' 



0. (57) 



where Greek indices are used for compound indices gath- 
ering cluster site, spin and possible band indices. 



are treated exactly. In fact, the Hamiltonian H is not 
altered in any way; the Weiss fields are introduced to let 
the variational principle act on a space of self-energies 
that includes the possibility of specific long-range orders, 
without imposing those orders. 

Steps towards a VCA calculation are as follows: 

1. Choose the Weiss fields to add, aided by intuition 
about the possible broken symmetries to expect. 

2. Set up a procedure to calculate the functional (55). 

3. Set up a procedure to optimize the functional, i.e., 
to find its saddle points, in the space of variational 
parameters. 

4. Calculate the properties of the model the saddle point. 



A. Practical calculation of the PotthofF functional 

Let x denote the (finite) set of variational parameters 
to be used. The Potthoff functional becomes the function 



TL 

n t (x) = n'-— EE lndet 



l-V(k)G'(k,w) (58) 



Once the cluster Green function is known by the meth- 
ods described in Sect. Ill, calculating the functional (58) 
requires an integral over frequencies and wavevectors of 
an expression that requires a few linear-algebraic opera- 
tions to evaluate. Two different methods have been used 
to compute these sums, described in what follows. We 
will see that the second method, entirely numerical, is 
much faster than the first one, which is partly analytic, 
a result that may seem paradoxical. 



1. Method I : Exact frequency integration 



V. THE VARIATIONAL CLUSTER 
APPROXIMATION 

The Variational Cluster Approximation 1 ' 22 (VCA), 
also called Variational Cluster Perturbation Theory 
(VCPT), can be viewed as an extension of Cluster Pertur- 
bation Theory in which some parameters of the cluster 
Hamiltonian are set according to Potthoff's variational 
principle through a search for saddle points of the func- 
tional (55). The cluster Hamiltonian H' is typically aug- 
mented by Weiss fields, such as the Neel field (44) that 
allow for broken symmetries that would otherwise be im- 
possible within a finite cluster. The hopping terms and 
chemical potential within H' may also be treated like ad- 
ditional variational parameters. In contrast with Mean- 
Field theory, these Weiss fields are not mean fields, in 
the sense that they do not coincide with the correspond- 
ing order parameters. The interaction part of H (or H') 
is not factorized in any way and short-range correlations 



The integral over frequencies in (55) may be done an- 
alytically, with the result 2 ' 5 

n(x) = n'(x)- xx + ^E E ^(k) (59) 

W r<° k W r (k)<0 

where the uj' r are the poles of the Green function G' 
in the Lehmann representation (B4) and the w r (k) are 
the poles of the VCA Green function (G " 1 (k) — Z) _1 . 
The latter are the eigenvalues of the R x R matrix 
L(k) = A + Q f V(k)Q (see Appendix B 1). R is the num- 
ber of columns of the Lehmann representation matrix Q, 
basically the total number of iterations performed in the 
band Lanczos procedure. 

In practice, the first sum in (59) is readily calculated. 
The second sum demands an integration over wavevec- 
tors. For each wavevector k, one must calculate L(k) and 
find its eigenvalues, a process of order R 3 . Other linear- 
algebraic manipulations leading to the diagonalization of 
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L(k) are typically less time-consuming than the diagonal- 
ization itself. The computation time therefore goes like 
N^R 3 , where is the number of points in a mesh cover- 
ing the reduced Brillouin zone (in fact half of the reduced 
Brillouin zone, since inversion symmetry is assumed). 



2. Method II : Numerical frequency integration 

An alternate method of computing the sums in (58) is 
to perform them in the reverse order, i.e., to first compute 
the wavevector sum for a fixed frequency u>, and then in- 
tegrating over frequencies numerically. The method used 
to sum over wavevectors is exactly the same as in Method 
I above : a wavevector mesh is set up in the reduced Bril- 
louin zone. This mesh is either a fixed, regular grid, or 
an adaptive mesh that is refined recursively as needed by 
comparing a two- and three-points Gauss-Legendre eval- 
uations within each cell (more accurately, the number of 
function evaluations in each cell is 2 d and 3 d , d being the 
dimension of space) . 

In the limit of zero temperature, the second term of 
the Potthoff functional (58) may be written as 



dco L 



^lndct [l - V(k)G'H 



(60) 



where the frequency integral / is carried along a closed, 
counterclockwise contour C that encloses the negative 
real axis, following the usual prescription with Green 
functions. We show in Appendix E that this integral 
reduces to 

PQO J T 

I = I VjV?H det(1 " V(£)G/(iic)) -^-rt 

(61) 

We refer to Appendix E for details, and for a discussion 
of the merits of this approach compared to the exact 
method above. This is the method that we generally 
follow. 
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FIG. 5: Potthoff functional as a function of Neel Weiss field 
M for various values of U, at half-filling, calculated on a 2 x 2 
cluster. The positions of the minima are indicated. 




FIG. 6: Optimal Neel Weiss field M and corresponding order 
parameter, as a function of U , at half-filling, calculated on 
a 2 x 2 cluster. Also shown is the ordering energy, i.e., the 
difference between the energy density of the normal state and 
that of the Neel state (in fact the difference between the grand 
potentials of the two solutions, since they both sit at half- 
filling). 



B. Example : Antiferromagnetism 

Let us start our examples with Neel antiferromag- 
netism. The corresponding Weiss field is defined in (44). 
Fig. 5 shows the Potthoff functional as a function of Neel 
Weiss field M for various values of U, at half-filling, 
calculated on a 2 x 2 cluster. We note three solutions 
per curve: two equivalent minima located symmetrically 
about M — 0, and a maximum at M = corresponding 
to the normal state solution. The normal and AF solu- 
tions both correspond to half-filling, and the AF solution 
has a lower energy density S = fi + fin. We therefore 
conclude, on this basis, that the system has AF long- 
range order. Note that, as U is increased, the profile of 
the curve is shallower and the minimum closer to zero. 
Indeed, for large U, the half-filled Hubbard model is well 



approximated by the Heisenberg model with exchange 
J = At 2 /U, and the curve should (and will) scale to- 
wards a fixed shape when Q/J is plotted against M/J 
(both dimensionless quantities). Fig. 6 shows how the 
optimal Weiss field and the Neel order parameter vary as 
a function of U. The Weiss field vanishes both as U — > 0, 
where the order disappears, and as U — * oo. In both 
limits the energy difference between normal and broken 
symmetry state (or 'condensation energy') goes to zero 
(Fig. 6), and so should the critical (Neel) temperature. 
The order parameter (Cm) increases monotonically with 
U and saturates. 

Fig. 7 shows the Potthoff functional as a function of 
Neel Weiss field M for various cluster sizes, at half- filling 
and U = 8. There is a clear and monotonous size depen- 
dence of the position of the minimum. In particular, the 
optimal Weiss field decreases as cluster size increases. 
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FIG. 7: Potthoff functional as a function of Neel Weiss field 
M for various cluster sizes, at half- filling and U — 8. The 
clusters used (from top to bottom) are: 2 x 2, 2 x 3, 2 x 4, 
B10 - see Fig. 1 -, and 3x4. The positions of the minima 
are indicated. 




0.6 

scaling factor 

FIG. 8: Top: Optimal Neel Weiss field for the half-filled Hub- 
bard model at U — 16, as a function of scaling parameter. 
Blue points : the scaling parameter is 1 — 1/L, and the scal- 
ing is poor. Red points : the scaling parameters is the number 
of cluster links divided by 2L - this takes open boundary con- 
ditions into account. We see how the Weiss field goes to zero 
in the thermodynamic limit. Bottom: Same, for the Neel 
order parameter, which tends to a finite value in the thermo- 
dynamic limit. Against the second scaling parameter works 
better. 



This should not worry us, quite on the contrary. The 
Weiss field is needed only because spontaneously broken 
symmetries cannot arise on a finite cluster. The bigger 
the cluster, the easier it is to break the symmetry and 
the optimal Weiss field should tends towards zero as the 
cluster size goes to infinity. Finite-size scaling is gen- 
erally very difficult, because cluster sizes are small and 
clusters vary in shape as well as size. Moreover, open 
boundary conditions are used rather than periodic ones, 
which adds edge effects to size effects. One needs to de- 
fine a scaling parameter q, ranging between and 1, that 
somehow defines the "quality" of the cluster (q = 1 being 
the thermodynamic limit). Fig. 8 shows the optimal Neel 
Weiss field as a function of two possibilities for the scal- 
ing factor q, for the half-filled Hubbard model at U = 16. 
The first possibility (blue dots) is q = 1 — 1/L, which 
does not take into account the shape of the cluster. The 
second possibility (red dots) corresponds to q defined as 
the number of links on the cluster, divided by twice the 
number of sites. This also goes to 1 in the thermody- 
namic limit (for the square lattice), but this time takes 
into account the boundary of the cluster. Indeed, 1 — q 
corresponds to the fraction of links of the lattice that are 
"inter-cluster" and thus treated "perturbatively" in the 
CPT sense. In that case, the scaling is good, as the op- 
timal Weiss fields extrapolates very close to zero in the 
q —> 1 limit. At the same time, the AF order parameter 
also decreases, bu extrapolates to a finite value, as shown 
on the same figure 



C. Superconductivity 

Superconductivity requires the use of pairing fields as 
Weiss fields, i.e., of operators creating Cooper pairs at 
specific locations. Generally, pairing fields have the form 



O sc = AijCiiCji + H.c 



(62) 



Different types of superconductivity correspond to differ- 
ent pairing functions A^-. For instance, ordinary (local) 
s-wave pairing (a la BCS) corresponds to = dij. On 
a square lattice, what is usually known as d x 2_ y 2 pairing 
corresponds to 



1 

-1 



if 
if 



- v, = ±x 



- tj = ±y 



whereas d xy pairing corresponds to 



1 if r,- Tj = ±(x + y) 
-1 if r, - v 3 = ±(x - y) 



(63) 



(64) 



The above two pairing are spin singlets. 

Pairing fields, once introduced in the cluster Hamilto- 
nian H' as Weiss fields, do not conserve particle number 
(but conserve spin). This increases the computational 
burden, since now the Hilbert space must be increased to 
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FIG. 9: Profile of the Potthoff functional as a function of 
Weiss field for various superconducting pairing fields. The 
extended s-wave is defined as the same as (63), but without 
the sign change between x and y directions. 



include all sectors of a given total spin. In practice, one 
uses the Nambu formalism, which in this case amounts 
to a particle-hole transformation for spin-down operators. 
Indeed, if we introduce the operators 

c a = c„t and d a = (65) 

then the pairing fields look like simple hopping terms 
between c and d electrons, and the whole cluster Hamil- 
tonian can be kept in the standard form (1), albeit with 
hybridization between c and d orbitals. 

Fig. 9 illustrates the dependence of the Potthoff func- 
tional on various superconducting pairing fields (gencr- 
ically denoted A). In that case, only d x 2_ y 2 pairing 
leads to a nontrivial solution. Others are piece-wise 
monotonously increasing or decreasing function, with a 
single zero-derivative point at A = 0. 

D. Thermodynamic consistency 

One of the main difficulties associated with VCA (or 
CPT) is the limited control over electron density. In the 
absence of pairing fields, electron number is conserved 
and clusters have a well-defined number of electrons. 
This makes a continuously varying electron density a bit 
hard to represent. Of course, one may simply vary the 
chemical potential fi and look at the corresponding vari- 
ation of the electron density, given by the trace of the 
Green function (schematically, TrG, see Appendix D). 
This provides a continuously varying estimate of the den- 
sity as a function of fi. An alternate way of estimating 
the density is to use the relation 
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FIG. 10: Comparisons of the estimates of the electron density 
n as a function of chemical potential fi, with different methods 
of calculation, for the normal solution, at U = 8, on a 2 x 2 
cluster. The subscript 'cons.' means that the corresponding 
quantities were computing in a thermodynamically consistent 
way, by using fi' as a variational parameter. 



where the grand potential £1 is approximated by the Pot- 
thoff functional at the solution found, and fi is varied 
as an external parameter. The problem is that the two 
estimates do not coincide (see Fig. 10). In other words, 
the approach is not thermodynamically consistent. The 
recipe to make it consistent is simple : the chemical po- 
tential fi' of the cluster should not be assumed to be 
the same as that of the lattice system (fi), but should 
be treated as a variational parameter. If this is done, 
then the two methods for calculating n given precisely the 
same result (see Fig. 10), and this can easily be proven 
in general. Results on a Hubbard model for the cuprates 
with thermodynamic consistency are shown on Fig. 11; 
see also Ref. 23. 



E. Searching for stationary points 

Let Xi be the n different variational parameters used 
in VCA. Once the function fi(x) may be efficiently cal- 
culated, it remains to find a stationary point of that 
function. This point is not necessarily a minimum in 
all directions. Indeed, experience has shown that to is 
a maximum as a function of the cluster chemical poten- 
tial fi' , while it is generally a minimum as a function of 
symmetry-breaking Weiss fields like M or A. 

The Newton-Raphson algorithm allows one to find sta- 
tionary points with a small number of function evalua- 
tions. One starts with a trial point x and an initial 
step h. Let denote the unit vector in the direction 
of axis i of the variational space. The function lo is 
then calculated at as many points as necessary to fit a 
quadractic form in the neighborhood of x . This requires 
(n+l)(n+2)/2 evaluations, at points like x , x ±hei, and 
a few of Xo + h(ei + e^). The stationary point Xi of that 
quadratic form is then used as a new starting point, the 
step h is reduced to a fraction of the difference |x x — x |, 
and the process is iterated until convergence on |x, — 1 1 
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FIG. 11: Order parameters for d x 2_ y 2 pairing and Neel an- 
tiferromagnetism for a model of the high-T c cuprates with 
U = 8, diagonal hopping ii = —0.3 and third neighbor hop- 
ping t2 — 0.2. Calculations are performed on a 3 x 4 cluster. 
Three solutions are displayed: (1) a pure d x 2_ y 2, obtained 
with two variational parameters (// and A x 2_ y 2); (2) a pure 
Neel solution obtained by varying and the Neel Weiss field 
M; a homogeneous coexistence solution obtained by varying 
fi' , M and A x 2_ y 2 (From M. Guillot, MSc thesis, Universite 
de Sherbrooke). 
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FIG. 12: Examples of clusters with baths. Bath sites are 
square, cluster sites blue circles. Bath parameters for the 
normal solution are indicated on the two top clusters, while 
conventional labels of orbitals are indicated for the (4+8)-site 
cluster. 



is achieved. A variant of this method, the quasi-Newton 
algorithm, may also be used, in which the full Hessian 
matrix of second derivatives is not calculated. It requires 
in general more iterations, but fewer function evaluations 
at each step. 

The advantage of the Newton-Raphson method lies in 
its economy of function evaluations, which are very ex- 
pensive here: each requires the solution of the cluster 
Hamiltonian. Its disadvantage is a lack of robustness. 
One has to be relatively close to the solution in order to 
converge towards it. But one typically runs parametric 
studies in which an external (i.e. non variational) param- 
eter of the model is varied, such as the chemical potential 
fi or the interaction strength U. In this context, the so- 
lution associated with the current value of the external 
parameter may be used as the starting point for the next 
value, and in this fashion, by proximity, one may conduct 
rather robust calculations. 

A more robust method, albeit more time consuming, 
is the conjugate-gradient algorithm, which we will not 
explain here as it is amply documented and fairly com- 
mon. However, this algorithm finds minima (or maxima), 
not saddle points in general. We must therefore take the 
extrinsic step of identifying parameters (like // above) 
that are expected to drive maxima of u>, and a comple- 
mentary set of parameters (like M and A above) that 
drive minima of to. One then, iteratively, finds maxima 
and minima with the two sets of parameters in succes- 
sion, and stops when convergence on |xj — Xj_i| has been 
achieved. This method is suitable to find a first solu- 
tion when the Newton-Raphson method fails to deliver 



one. It may however converge to minima that are in 
fact singularities of lo, i.e., points where the derivatives 
are not defined. Such points may occur as the result of 
energy-level crossings in clusters and are an artifact of 
the finite-cluster size. 



VI. THE CELLULAR DYNAMICAL MEAN 
FIELD THEORY 

The Cellular dynamical mean-field theory (CDMFT) 
- also called Cluster dynamical mean-field theory - 
is a cluster extension of Dynamical mean-field theory 
(DMFT). Since there is no real pedagogical gain in de- 
scribing first DMFT, we will proceed directly to CDMFT, 
in the context of a an exact diagonalization solver. 

The basic idea behind CDMFT is to model the effect 
on the cluster of the remaining degrees of freedom of the 
lattice by a bath of uncorrelatcd orbitals that exchange 
electrons with the cluster, and whose parameters are set 
in a self-consistent way. Explicitly, the cluster Hamilto- 
nian H' takes the form 

H' = - ^2 t^c^c v + U ^2 n ^n ai 

(67) 

where a a annihilates an electron on a bath orbital la- 
belled a. The label a includes both an 'bath site' index 
and a spin index for that 'site'. The bath is characterized 
by the energy of each orbital (e Q ) and by the bath-cluster 
hybridization matrix 6^ a (the index /i includes cluster 
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site, spin and band indices). This representation of the 
environment through an Anderson impurity model was 
introduced in Ref. 24 in the context of DMFT (i.e., a 
single-site cluster). Note that 'bath site' is a misnomer, 
as bath orbitals have no position assigned to them. 

The effect of the bath on the electron Green function 
is encapsulated in the so-called hybridization function 



rV(w) 



^ W - £„ 



(68) 



which enters the electron Green function as 

G' _1 =w-t- V(oj) - E(w) (69) 

This is shown in Appendix F in the non-interacting case 
(Y. = 0). By definition, the only effect of adding the 
electron-electron interaction is to add the self-energy Z, 
as above. 

A. Bath degrees of freedom and SFA 



0.95 




0.85 



FIG. 13: Electron density as a function of chemical potential, 
for the one-dimensional Hubbard model, at U = 3. The black 
curve is the exact result from the Lieb-Wu solution using the 
Bethe Ansatz. The red curve is the SFA calculation on the 
2-site cluster with 4 bath sites shown on Fig. 12. The other 
curves are obtained with the CDMFT algorithm (parameters 
explained in the text). 



The CDMFT Hamiltonian (67) defines a valid ref- 
erence system for Potthoff's self-energy functional ap- 
proach, since it shares the same interaction part as the 
lattice Hamiltonian H and since each cluster of the su- 
perlattice has its own identical, independent copy. From 
the SFA point of view, the bath parameters {e a , 9fj, a } can 
in principle be chosen in such a way as to make the Pot- 
thoff functional stationary. A subtlety arises: the bath 
system must be considered part of the original Hamilto- 
nian H, albeit without hybridization to the cluster sites, 
in order for both Hamiltonians to describe the same de- 
grees of freedom; but within H we are free to give the 
bath trivial parameters (e a = 0). Performing VCA-like 
calculations with bath degrees of freedom is illustrated 
in Ref. 25, and on Fig. (13) below. 

When evaluating the Potthoff functional in the pres- 
ence of a bath, one must add a contribution from the 
bath to Tr ln(— G'), which takes the form 
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FIG. 14: Two of the four bath parameters used in the calcu- 
lations of Fig. 13, as a function of /i. The legend is the same 
as for Fig. 13. 



^bath = ^2 £a ( 70 ) 

and which comes from the zeros of the cluster Green func- 
tion induced by the poles of the hybridization function. 
Note that the zeros coming from the self-energy cancel 
out in Eq. (59) between the contribution of Tr ln(— G') 
and that of Tr ln(— G), but not those coming from T(w), 
as they only occur in G'. 



B. The CDMFT self-consistent procedure 

However, in practice, CDMFT does not proceed in this 
way, i.e., it does not look for a strict solution of the Euler 



equation (57) . It tries instead to set each of the terms be- 
tween brackets to zero separately. Since the Euler equa- 
tion (57) can be seen as a scalar product, CDMFT re- 
quires that the modulus of one of the vectors vanish to 
make the scalar product vanish. From a heuristic point 
of view, it is as if each component of the Green function 
in the cluster were equal to the corresponding compo- 
nent deduced from the lattice Green function. Clearly, 
the left-hand side of Eq. (57) cannot vanish separately 
for each frequency, since the number of degrees of free- 
dom in the bath is insufficient. Instead, one adopts the 
following self-consistent scheme (see Fig. 15): 

1. Start with a guess value of the bath parameters 
(0^a,£a), that define the hybridization function (68). 
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of the order of 20 or 40 (in units of i" 1 ). Even when the 
total number of cluster plus bath sites in CDMFT equals 
the number of sites in a VCA calculation, CDMFT is 
much faster than the VCA since the minimization of a 
grand potential functional requires many exact diagonal- 
izations of the cluster Hamiltonian H' . 
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FIG. 15: The CDMFT algorithm with an exact diagonaliza- 
tion solver. 



2. Calculate the cluster Green function G(lu) with the 
Exact diagonalization solver. 



3. Calculate the super lattice-averaged Green function 

(71) 



6(0,) = ^-^ 



and the combination 

%"V) = G- X + E(u;) (72) 

4. Minimize the following distance function: 

d= l^ + M-f-r^-Sfo- 1 )^! 2 (73) 

over the set of bath parameters (changing the bath 
parameters at this step does not require a new solution 
of the Hamiltonian H' , but merely a recalculation of 
the hybridization function T). 

5. Go back to step (2) with the new bath parameters 
obtained from this minimization, until they are con- 
verged. 

In practice, the distance function (73) can take 
various forms, for instance by adding a frequency- 
dependent weight in order to emphasize low-frequency 
properties 26-28 or by using a sharp frequency cutoff. 29 
These weighting factors can be considered as rough ap- 
proximations for the missing factor 5H'(u>)/5t' in the 
Euler equation (57). The frequencies are summed over 
on a discrete, regular grid along the imaginary axis, de- 
fined by some fictitious inverse temperature (3, typically 



C. Examples 

Let us start with a one-dimensional example. Fig. (13) 
illustrates the variability of CDMFT results related to 
the choice of the distance function. The cluster used has 
two sites and four bath sites (see Fig. 12), and the vari- 
ous curves represent the electron density n as a function 
of chemical potential /i for the one-dimensional Hubbard 
model at U — 3. The exact results, from the Lieb-Wu 
solution, is shown in black, as well as the SFA result 
coming from an exact solution of the Euler equation (57) 
for that system. The four CDMFT results shown differ 
by the value of (3 and that of the sharp frequency cut- 
off lo c . In addition, one of the curves was obtained by 
weighing the different frequencies by a factor An 
important characteristic of the exact result is that infi- 
nite compressibility dn/dfj, at the point where the gap 
ends, i.e., when the density curve hits n = 1, at a value 
fj, c (U) of the chemical potential. The SFA and CDMFT 
do quite well in accounting for the infinite compressibil- 
ity, contrary to other approaches (e.g. one-site DMFT). 
On this small two-site system, they do not find the cor- 
rect (U c , but increasing the bath size would improve on 
this. By playing with the distance function, on may bring 
the curves closer to or further from the exact result, but 
there is no guarantee that the most successful distance 
function in this case will be as profitable when the exact 
solution is unknown! In principle, the SFA curve (red) is 
the one that best represents what can be achieved with 
this system, and the various CDMFT curves are to be 
judged against not the exact result, but against the SFA 
curve. 

Next, consider the two-dimensional cluster illustrated 
in the lower part of Fig. 12. This 4-site, 8-bath site cluster 
is the main cluster used in CDMFT simulations of high- 
ly cuprates using the two-dimensional Hubbard model. 
It is useful in that case to view the orbitals numbered 
5 to 8 as a first bath set, and the orbitals numbered 9 
to 12 as a second bath. Each site of the cluster is con- 
nected to one orbital of each set. In studying the normal 
state, and taking into account the symmetries of the clus- 
ter, we would need 4 bath parameters: one bath-cluster 
hopping and one bath energy for each set. In order to 
treat a possible antiferromagnetic phase, one must mod- 
ify the bath energies and hopping in a spin-dependent 
way. The grey and white squares on the figure then dis- 
tinguish orbitals of a given bath according to their shift 
in site energy (of opposite signs for opposite spins). The 
corresponding bath-cluster hybridization may also be dif- 
ferent, which makes a total of 8 parameters. Finally, in 
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FIG. 16: Neel (AF) and d-wave (dSC) order parameters ob- 
tained from CDMFT applied to the (4+8)-cluster of Fig. 12, 
for the two-dimensional Hubbard model with U — 8 and di- 
agonal hopping t' = —0.3. The data is shown as a function 
of the calculated lattice density n. The order parameters are 
calculated using the same operators as in the corresponding 
VCA calculation illustrated on Fig. 9, even though these op- 
erators played no role in the solution: they are merely used 
as a probe. In this calculation, we set (3 — 20 and a sharp 
cutoff Lu c = 3 was used. The dSC and AF solutions were 
both allowed simultaneously (9 bath parameters) and there 
are regions of coexistence of the two orders. 



order to study <i-wave superconductivity, we introduce 
pairing within each bath (red dotted lines on the figure) , 
vertical and horizontal pairing being of opposite signs. 
This introduces an additional parameter, for a total of 
9. At this point, an important remark is in order : For- 
mula (68) for the hybridization function only applies if 
the bath orbitals are not hybridized between themselves. 
The (i-wave pairing just described certainly breaks that 
condition. This is not a problem, however, if we perform 
a change of variables within bath degrees of freedom (a 
Bogoliubov transformation) prior to solving the problem 
numerically, such as to make the bath Hamiltonian di- 
agonal. Then the poles of the hybridization function no 
longer correspond to the bath energies as defined origi- 
nally in the model, but rather to the eigenvalues of the 
bath Hamiltonian. 

Results of a CDMFT calculation on this system are 
shown in Fig. 16. Comparing with the VCA result of 
Fig. 11, we notice first the similarities: the existence of 
a dSC phase away from half-filling for both electron and 
hole doping and the possibility of homogeneous coexis- 
tence between antiferromagnetism and d-wave supercon- 
ductivity. But differences are obvious : the VCA diagram 
is more asymmetric than the CDMFT one in terms of 
electron vs hole doping. Both calculations agree on the 
critical doping for antiferromagnetism on the hole-doped 
side (~ 10%), but not on the electron-doped side. The 
VCA result does not show homogeneous coexistence be- 
tween AF and dSC on the hole-doped side - although it 
appears on smaller clusters. At this point it is not clear 
whether these differences arise because of the methods 



themselves rather that the particular way they were ap- 
plied (choice of Weiss fields, bath configuration, distance 
function, etc.). In particular, the exact SFA result for 
the system used in CDMFT has not yet been calculated. 



APPENDIX A: CLUSTERS AND KINEMATICS 

In this appendix we will review the kinematics of 
cluster decompositions, and introduce notation used 
throughout this paper. The spatial dimension D of the 
lattice will be left general. 

Cluster methods are based on a cluster decompostion 
of the model, i.e., on a tiling of the original lattice 7 
with identical clusters of L sites each. Mathematically, 
this corresponds to introducing a superlattice T, whose 
sites form a subset of the lattice 7 and will be labelled 
by latin indices such as (m, n...). This superlattice is 
generated by D basis vectors ...d belonging to 7, i.e., 
every site r m of the superlattice may be expressed as 
an integer linear combination of these basis vectors. As- 
sociated with each site of the T is a cluster of L sites, 
whose shape is not uniquely determined by the superlat- 
tice structure. The sites of the clusters will be labelled by 
latin indices (a, 6, . . . ). Each site of the original lattice 
7 can be expressed in a unique way as a combination of a 
superlattice vector r m and of a site r a within the cluster: 
i"; = r m +r a . This means that the index i may in fact be 
replaced by a compound index (m, a), and we have the 
following equivalence between summations: 
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The number of sites in the cluster is simply the ratio of 
the unit cells volumes of the two lattices. In D = 3, this 
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(the above formulae can be adapted to D = 2 by setting 
e 3 = (0,0,1)). 

The Brillouin zone of the original lattice, denoted BZ 7 , 
contains L points belonging to the reciprocal superlattice 
r*. The Brillouin zone of the superlattice, BZp, has a 
volume L times smaller than that of the original Brillouin 
zone. Any wavevector k of the original Brillouin zone can 
be uniquely expressed as 



k = K + k 



(A3) 



where K belongs both to the reciprocal superlattice and 
to BZ 7 , and k belongs to BZr (see Fig. 17). Thus, we 
have the equivalent summations 
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The passage between momentum space and real space, 
by discrete Fourier transforms, can be done either di- 
rectly (i <-> k), or independently for cluster and super- 
lattice indices (m <-> k and a <-> Q): 

Si = E e ^/( k ) /( k ) = ^E e-<k,p Vj 

/ m = £V fcr -/(k) /(k) = ^$>- <kr »/ m 

k m 



K 



(A5) 

where / stands for a generic one-index quantity. Quasi 
continuous indices, like k and k, are most of the time 
indicated between parentheses. This notation may right- 
fully be deemed capricious, since the labels i and k take 
the same number N of values, but we adopt it nonethe- 
less as it helps reminding us that the values of the labels 
are closely separated. These discrete Fourier transforms 
close by virtue of the following identities 
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where 5j is the usual Kronecker delta, used for all labels 
(since they are all discrete): 

1 if a = 



S a = 



S a f3 = S a -/3 > (A9) 



otherwise 
and the A's are the so-called Laue functions: 



A 7 (k) = <5 k+Q 
Ar(k) = J2 S * 



'k+P 



(A10) 
(All) 
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Laue functions are used instead of Kronecker deltas in 
momentum space because of the possibility of Umklapp 
processes. Note especially that even though 

S k = S k S K (k = k + K), (A12) 

the same does not hold for the Laue functions: 

A 7 (k) ^ A r (k)A 7 (K) . (A13) 

Instead we have the following relations: 

A r (k) =J]A 7 (k + K) (A14) 

K 

A 7 (k) = A 7 (k + K) = <%A 7 (K) (A15) 




(Tj 71") 



FIG. 17: The reduced Brillouin BZr zone associated with the 
10-site cluster of Fig. 1. A wavevector k has a unique decom- 
position k = k + Q, where Q is one of the L elements of the 
reciprocal superlattice that belongs to the original Brillouin 
zone BZ-v. 



which reflect the arbitrariness in the choice of Brillouin 
zone of the superlattice (we use the term Brillouin zone 
in a rather liberal manner, as a complete and irreducible 
set of wavevectors, and not as the Wigner-Seitz cell of 
the reciprocal lattice.) 

A one-index quantity like the destruction operator Cj — 
c m ,a can be represented in a variety of ways, through 
partial Fourier transforms: 

L ^ ' (A16) 



c ^ k ) = A7E e " k ' rmc " 



c m.K 



~ TV ^ ,a 



c K (k) = 

c ( k )-lrE e " lk ' ri ^ 



(A17) 
(A18) 
(A19) 



The last two representations are not identical, since the 
phases in the two cases differ by k • r a . They are related 
by a L-dimensional unitary matrix: 



-(k + K) = ]TAKK'(k)cK(k) 



where 



K' 



A K K<(k)^£ e ~ ir " (k+K ~ K ' 



(A20) 



(A21) 



A two-index quantity like the hopping matrix iy may 
thus have a number of different representations. Due 
to translation invariance on the lattice, this matrix is 
diagonal when expressed in momentum space: t(k, k') — 
e(k)<5k,k', £k being the dispersion relation: 



k 



(A22) 
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However, we will very often use the mixed representation 

i = (0,a) 



j = (m, b) 



(A23) 



For instance, if we tile the one-dimensional lattice with 
clusters of length L = 2, the nearest-neighbor hopping 
matrix, corresponding to the dispersion relation e(k) = 
— 2tcos(fc), has the following mixed representation: 



t(fc) 




-2ik 



(A24) 



APPENDIX B: LEHMANN REPRESENTATION 
OF THE GREEN FUNCTION 

By inserting a completeness relations in the expression 
(26) for the zero-temperature Green function, one finds 
the Lehmann representation: 



m 



1 



(m|4|fi) 



OJ — E n 

1 



En 



oj + E n — Eq 



(n|c p |Q) 



(Bl) 



The two sums are over different sets of eigenstates, in the 
spaces with one more and one less electron, respectively. 
Let us introduce the notation 



QW = <fi|c„|m> 



Q$ = <fi|cj» (B2) 



as well as uiff — E m — E > and oj^ = -E n + E < 
to write 



<-*»Au) - 2^ fir + 2^ 



(e) 1 (ft) 
Uf — OJ m «. W — W„ 



(B3) 



The Q)fm form a 2L x 7V< e ) matrix, where A^( e ) is the 
number of states \m) that give a nonzero contribution to 

the first sum above. Likewise, The Q^X form a 2L x TV^ 1 ) 
matrix. Let ./V = A^ e ) + and let us introduce a 

2L x iV matrix Q by joining vertically the matrix 
below the matrix Q^ e \ and let w r denote the elements of 
the concatenated sets {w™} and {oj^}. Then we can 
write 



^ OJ — 0J r 



(B4) 



If we introduce the diagonal matrix A rs = <5 rs o; r and 

1 



- A 

then we have the matrix expression 
G(w) - Qg(^)Q t 



(B5) 



(B6) 



This is a very general representation of the exact cluster 
Green function. 



1. The Lehmann representation and the CPT 
Green function 



Let us see how the CPT Green function can be explic- 
itly represented in terms of the Lehmann representation 
(B4). The CPT Green function (3) can be written as 23 



G(k,w) = 



(QgHQt)-i-V(k) 
= QgMQt + (QgHQt)V(QgHQ+) + 

= Q(gM + g(w)(QtVQ)g(w) + ---)Qt 



= Q 



l 



oj - L(k) 



(B7) 



where L(k) = A + QtV(k)Q. The poles of G(k,w) are 
those of [oj — L(k)] -1 , which we denote as w r (k). They 
are simply the eigenvalue of the N x N matrix L(k). 
Let U(k) the matrix that diagonalizes L(k), such that 



U(k)L(k)U f (k) = A(k) 
where A(k) is diagonal. Then we write 



(B8) 



G(k,c) = Q * Qt = QU(k) l 7f -(QU(k))t 

oj - L(k) oj — A(k) 

(B9) 

which again is of the same form as (B4) , with Q replaced 
by Q(k) = QU(k). 

The representations (B4) or (B9) ensure the positivity 
of the cluster Green function and the CPT Green func- 
tion respectively, i.e., the positive character of the corre- 
sponding spectral functions. Indeed, the local (cluster) 
spectral weight is 



A a (oj) = -2 lim Llm G' aa (oj + i V ) 

7?— >0 



and 



\Qa 



(BIO) 



(Bll) 



This expression has poles on the real axis only with pos- 
itive residues, and this garantees that the corresponding 
spectral function A a (oj) is positive. Moreover, the prop- 
erty QQ f = 1 ensures it is normalized. 

Positivity will also hold for the Fourier transformed 
(i.e. momentum) Green function. Indeed, let us consider 



G(k,,) = |^c- k '( r »- r ')G at 



a,f3 r 

£|W(k)|VH 



r br9r (uj) (B12) 
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where 



Vv(k) 



-ik-r 



(B13) 



The above expression is also manifestly positive. Of 
course we could have used the CPT Green function in the 
Fourier transform instead of the cluster Green function 
and the conclusion would be the same: it is a consequence 
of the Lehmann representation. 



APPENDIX C: GROUP THEORETICAL 
CONCEPTS 

This short appendix summarizes some key group- 
theoretical concepts necessary to understand the discus- 
sion of Sect. HIE. Of course this is no substitute to a 
text on group theory. It merely serves as a reminder to 
those who have some knowledge of it, or indicates some 
important concepts to those who don't. 

Let 25 denote the discrete symmetry group of the sys- 
tem and |® | the number of elements in the group (the 
order of the group). Elements of the group will be de- 
noted by latin letters like 3, h, etc. and gh will stand for 
group multiplication, i.e., the symmetry operation ob- 
tained by applying first h, then 3. Recall that the set © 
forms a group if the following conditions are met: 

1. The set must be closed under the group multiplication, 
i.e., if gi and 32 belong to ©, so must g\gi- 

2. There must be a neutral element e (the identity trans- 
formation) such that eg = ge. 

3. Each element g must have a unique inverse g~ x such 
that gg- 1 =g- 1 g = e. 



4. The group operation must be associative : (31(72)33 = 
5i(5233)- 

The group multiplication may or may not be commuta- 
tive. In the first case, the group is said to be Abelian. 

The simplest non trivial group is C2 , the group of two 
elements formed by the identity transformation and a 
7r rotation (or, equivalently, an inversion). Examples of 
cluster systems with this symmetry group are illustrated 
on Fig. 18. Another common symmetry group is C^v, 
which consists of a 7r-rotation ci and two unequivalcnt 
reflexions (p\ and (72)- This is the symmetry group of 
a rectangular cluster, or of a square cluster with d x 2_ y 2 
pairing, for instance. Examples are illustrated on Fig. 19. 

A representation of the group is a set of matrices 
that behave exactly like the group elements when group 
multiplication is mapped onto matrix multiplication (i.e., 
there is an isomorphism between the abstract group and 
the set of matrices). In practice, quantum mechanics 
deals with group representations. The word representa- 
tion is also applied to the vector space (or module) on 
which the matrix representation is based. A representa- 
tion is said to be reducible if a change of basis can bring all 




FIG. 18: Clusters with C2 symmetry. Black and grey dots 
represent unequivalent sites, e.g., because of the presence of 
a Neel Weiss field. The top and bottom clusters are invari- 
ant under 7r-rotations, and the right-most cluster is invariant 
under a left-right inversion. The distinction is of course irrel- 
evant for the middle, one-dimensional cluster. 




FIG. 19: Clusters with Ci v symmetry. Dashed links on the 
square cluster illustrate the presence of a d x 2_ y 2 pairing field, 
which makes horizontal and vertical links unequivalent. 



group elements to the same block-diagonal form. Thus, 
reducible representations are direct sums of irreducible 
representations. It it the latter that are important, in 
great part because of Schur's lemmas, which imply that 
if a Hamiltonian matrix H commutes with all the group 
elements and if the basis states are arranged into irre- 
ducible representations, then H has no matrix elements 
between states belonging to different representations, i.e., 
it is block diagonal. We often say irrep for 'irreducible 
representation'. 

Two group elements 31 and 32 are said to be conjugate 
to each other if 31 = hr^gih for some element h of the 
group. This property is transitive, and therefore all the 
elements of a group may be organized into equivalence 
classes called conjugacy classes. Because elements of a 
conjugacy class are related by a similarity transforma- 
tion, they all have the same trace in a given represen- 
tation. This trace is called the character (denoted \) of 
the class in the said representation. The identity element 
e forms a conjugacy class all by itself, and its character 
is the dimension of the representation. It can be shown 
that the number of unequivalent irreps is the same as the 
number of conjugacy classes. Characters are often dis- 
played in tables (Fig. 20), as a function of the conjugacy 
class (horizontal) and irreps (vertical). These tables are 
extremely useful, for instance, to reduce tensor products 
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c 2 


e C2 


A 
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1 1 

1 -1 



C2v 


e 


C2 


<Tl 


0"2 


A! 


1 


1 
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1 


A 2 


1 


1 


-1 


-1 


Bi 


1 


-1 


1 


-1 


B 2 


1 


-1 


-1 


1 





e 


C2 


2c 4 


2<T1 


2(72 


Ai 


1 


1 


1 


1 


1 


A 2 


1 


1 


1 


-1 


-1 


Bi 


1 


1 


-1 


1 


-1 


B 2 


1 


1 


-1 


-1 


1 


E 


2 


-2 












FIG. 20: Character tables for the groups C2, Civ and Civ 



of irreps. Indeed, the trace of a matrix tensor product 
is the product of the traces of the factors, whereas the 
trace of a direct sum of matrices is the sum of the traces. 

Consider, for instance, the character table of C 2 „. We 
learn from it that this group has 4 disinct irreps, all of 
dimension one. Representation A\ is the trivial repre- 
sentation, with states even under all symmetry transfor- 
mations. A 2 contains states that are odd under cither 
reflexion. B\ and B2 contain states that are odd under 
a 7r-rotation, and under one of the two reflexions. C± v , 
on the other hand, has 5 irreps, of which the last (E) 
is two-dimensional and contains states that arc mixed or 
interchanged under 7r/2 rotations or reflexions. 

When working in a reducible space that is the direct 
sum of different irreducible representations, it is possible 
to project onto the various irreps making up that space 
using the following projection operators: 



p(a) _ 



E4 a) * 

g 



(CI) 



where g stands for a symmetry operation in the reducible 
space considered and \g 1S the character in the irrep a 
of the corresponding group element. 



APPENDIX D: CALCULATING AVERAGES 
FROM THE GREEN FUNCTION 

This appendix explains how to calculate the expecta- 
tion value of a one-body operator from the Green func- 
tion. 

A general one-body term of the Hamiltonian si written 

as 



O = SapC^Cp 



(Dl) 



where the indices a and (3 stand for all degrees of freedom: 
lattice site, spin and band. In the simple case of the 
number of electrons, the matrix s is diagonal: 



(D2) 



In the case of the antiferromagnetic order parameter, it 
is also diagonal and has the form 

Siaja' = %<W (-l) CT e^- r * Q = (tt, tt) (D3) 

We are interested in the expectation value density O — 

s a p{cicfi)/N. 

From the Lehmann representation of the Green func- 
tion, we see that (c^c^) is given by the integral of the 
Green function along a contour C< surrounding the neg- 
ative real frequency axis counterclockwise: 



(4 c «) = / 



dz 



-.G a p(uj) 



(D4) 



Therefore the expectation value we are looking for is 

dz 



O 



1 

N 



Sf3a (CgC a ) 



nJ c< 



2ni 



jtr [sGM 



(D5) 



Here the trace includes a sum over lattice sites, spin and 
band. In a mixed representation, with cluster sites in- 
dices and reduced wavevector instead of the lattice site 
index, this becomes 



O 



1 

N 



E 



c 2iri 



tr 



s(k)G(k,w) 



(D6) 



where we assumed that the matrix s is diagonal in k. 

Next, let us consider the asymptotic behavior of the 
Green function as z — > 00: G(w) — > \jz. This allows us 
to modify Eq. (D7) as follows: 



O 



TV ^ 

k 



f dz 






lc < 27ri 


tr 


s(k)G(k,w) 



trs(k) 



(D7) 

where p > (in practice, we use p ~ 1). The term 
we added does not contribute, since its unique pole lies 
outside of the contour. However, this term modifies the 
asymptotic behavior of the integrand, which nows decays 
as 1 j z 1 . This allows us to replace the contour C< by an 
integral along the imaginary axis, plus an infinite semi- 
circle that does not contribute, since the integrand falls 
faster than 1/z. 

Next, consider the part of the contour C< that lies 
above the real axis, and let us follow this contour clock- 
wise and call it C. Let C be the mirror image of C below 
the real axis, followed counterclockwise. To each z and 
dz of C correspond the mirror images z* and dz* on C, 
so that 

I[C] - / dz/M = [ dz*f(z*) (D8) 
Jc Jc 

If, in addition, the integrand is such that f(z*) = 
then 
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The integral of /(w) along the counterclockwise contour 
C< would then be 

J[C<] = I[C] - J[C] = 7*[C] - 7[C] = -2*Im 7[C] 

(D10) 

One of the properties of the Green function is its her- 
miticity: G a p(u)*) — G*^ a (uj). In the mixed Fourier 

representation, this is rather expressed as G(k, z*) — 
G^(— k, at). We also assume that s is Hermitian: s(k) = 
s^(— k) so that the expectation value is real. This means 
that the integrand of the expectation value respects the 
condition f(w*) — 

Finally, the expectation value has the expression 



O 



1 

N 



E 



duj 



Re i tr 



s(k)G(k,iw) 



trs(k) \ 
iuj — p J 
(Dll) 



APPENDIX E: EVALUATION OF THE 
POTTHOFF FUNCTIONAL 

In this appendix, we show how to evaluate the fre- 
quency integral (60). Let us take that contour to be the 
whole imaginary axis between —iR and iR, with a half- 
circle of radius R closing the contour on the left part of 
the complex plane. Let us first see what the behavior 
of the integrand is as \w\ — > oo. Since G'(w) ~ 1/oj as 
u — » oo, one may write, in that limit, 

lndet [l - V(k)G'(a;)l = tr ln \ l - V(k)G'H 



tr ln 



1 - 



V(k) 



trV(k) 



(El) 



The integration over wavevectors yields 



L £trV(k)=^5>-£C ( W ^oo) (E2) 



N 



The only contribution from these terms is the chemical 
potential, by definition. Thus the l/u> term in the inte- 
gral I is 



A(w) = 2L(ji-fi')- 



(E3) 



This term by itself does not contribute to the frequency 
integral. Indeed, let us arrange for the contour C, which 
normally should cross the origin, to avoid it along an 
infinitesimal semi-circle C\ of radius r\ centered at the 
origin and lying on the left-hand half- plane. This slight 
modification should not change the value of f2, if we re- 
fer to the exact method of the previous subsection, as 
the zero frequency does not contribute. Then the above 
term does not contribute, since the pole at u> = lies 
outside the contour. It can therefore be subtracted and 
the frequency integral has the following expression: 



L 




-5>det(l-V(k)G'M)-,4M 



in which the integrand now falls like 1/oj 2 at large fre- 
quencies, and the half-circle of radius R — > oo does not 
contribute to the integral. Let us use the properties 



G(— ix) = G(ix)* 



V(-ia;,k) = V(iz,-k)* (E4) 



to express the integral over the whole imaginary axis as 
an integral over the 'positive' imaginary axis only. With 
w = ix, we write 



dx L 

2~7T TV 



[lndct(l - V(k)G'(ix)) 



+ lndet(l- V(-k)*G'(ia;)*) 



Ci 



duj 
2wi 



(E5) 



Note that A(w) does not contribute to the integral along 
the imaginary axis, since it is odd in x; in other words, 
the principal value is taken and corresponds to the contri- 
bution of the small half-circle C\. We also note that the 
main part of the integrand does not have a contribution 
along the contour C\, since the logarithmic singularity 
near u — is integrable, i.e., leads to a vanishing con- 
tribution on Ci as 7] — > 0. Since integrating over k and 
over — k are equivalent, one further simplifies to 

r°° dr T I 

/ = I -^£ln|det(l-V(k)G'(iaO) 

(E6) 

The frequency integral is done by dividing the posi- 
tive imaginary frequency axis in three segments: [0,Ai], 
[Ai,A2] and [A2,oo). The constant Ai is a low-energy 
scale in the problem, like the lowest eigenvalue uj' r (up to 
some minimum), whereas A2 is a high-energy scale, like 
the largest eigenvalue ui' r (up to some maximum). On 
each segment a 20-point Gauss-Legendre integration is 
used. The last segment is in fact treated like an inte- 
gral over u = l/u> from to A^ 1 . The frequency inte- 
gral / dui f(oj) thus takes the form of a weighted sum 
J2 n Pnf{w n ), where f(u>) is the wavevector sum con- 
ducted at frequency ui. The required accuracy of the 
wavevector sum, which sets the number of points of the 
wavevector mesh, is conditioned by the weight p n (i.e. it 
does not have to be so large when p n is small). In ad- 
dition, the integrand of (58) may have sharp structures 
- thus requiring a fine wavevector mesh - at frequencies 
close to the real axis (the poles of G are all real) but is in- 
creasingly smooth as one moves away from the real axis. 
Thus the wavevector mesh may be redefined from time 
to time as one progresses along the imaginary frequency 
axis, making the mesh coarser and coarser without cost 
in accuracy. 

Tests have been conducted in order to compare the 
speed and accuracies of the two methods: the analytic 
frequency integration (AI) described in Sect. VA1, and 
the numerical frequency integration (NI) described in this 
section. In each fixed mesh and an adaptive mesh 

have been used. The results are displayed in Figs. 21 and 
22. 
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FIG. 21: Comparisons of the different integration methods in 
execution time and accuracy for a 4-site cluster. See text for 
details. 
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FIG. 22: Same as Fig. 21 for a 12-site cluster. 



Fig. 21 shows the relation between execution time (in 
seconds) on a 2.16 MHz Intel Core 2 Duo processor and 
the value of the functional obtained with four different 
methods: (i) Numerical integration (Nl) with a fixed 
wavevector grid, (ii) Nl with an adaptive wavevector grid, 
(iii) Analytic integration (Al) with a fixed wavevector 
grid and (iv) Al with an adaptive wavevector grid. For 
each method, different points correspond to different re- 
quired accuracies of the wavevector integral. The sys- 
tem used in that case was the two-dimensional Hubbard 
model with nearest-neighbor hopping t = 1, Coulomb 
repulsion U = 8 and chemical potential fi = 1, on a 
4-site cluster. The time to perform the exact diago- 
nalization is negligible here. The value of the chemi- 
cal potential used makes the system metallic, i.e., there 
are poles of the VCA Green function at zero frequen- 
cies, which is more demanding on the integral because 
of sharp features of the integrand. Deep within ordered 
phases there are generally gaps in the physical spectrum 
that make the integrals converge faster, but the location 
of phase boundaries, where these gaps disappear, is gen- 
erally of great interest. One sees from Fig. 21 that the 
Nl method with adaptive mesh converges fastest, about 
3 times faster than the Al method. The two horizontal 
red lines represent the two converged values (for Nl and 
Al). They differ by less than 10~ 6 , a difference due to the 
finite mesh used in the frequency integral. This accuracy 
is more than adequate for applications. 

Fig. 22 shows the same, this time for a 12-site clus- 
ter. The exact diagonalization is done beforehand, and 
is the same for all integration methods shown. In that 
case the number R of poles was 312, instead of R = 32 
for the 4-site cluster used in Fig. 21. Here the gain in 
using Nl is greater even, as the method is over 500 times 
faster than Al with a adaptive mesh. Indeed, it is pro- 
hibitively expensive to use the Al method for anything 
but very small clusters. Method I (Al) is of order N^R 3 , 
whereas Method II (Nl) involves linear-algebraic opera- 



tions on the Green function, of order NkN^L 3 (say, for a 
fixed wavevector grid and N u frequencies). The ratio of 
execution times between the two method should roughly 
be 

timc(AI) _ /R\ 3 1 

time(NI) ~ \LJ [ 1 

In the case of Fig. 21, R = 32, L = 4 and N u = 60, which 
gives a ratio of ~ 8. In the case of Fig. 22, R = 312, 
L = 12 and N u = 60, which gives a ratio of ~ 290. 
In both cases this overestimates by a factor 2 to 3 the 
measured advantage of Method II, which is nevertheless 
enormous. 



APPENDIX F: FERMIONIC BATHS AND 
HYBRIDIZATION FUNCTIONS 

In this short appendix we consider the Hamiltonian 

( F1 ) 

and show that the Green function obtained by tracing 
over the bath degrees of freedom has the form 

(G-V^-V-E^ 1 . (F2) 

a 

First of all, the full Green function associated with the 
above one-body Hamiltonian is 

GfunM - ^ (F3) 
where the full hopping matrix has the block form 
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where t is the hopping matrix within cluster degrees of 
freedom only, the hopping matrix between bath and 
cluster orbit als, and £ the diagonal matrix of bath ener- 
gies e a . The Green function obtained by tracing out the 
bath degrees of freedom is simply the restriction of Gf u n 
(and not of its inverse) to the cluster degrees of freedom 
only. The mathematical problem at hand is simply to 
invert a 2 x 2 block matrix 

i^A 2 i A22J (#21 £?22^ ' ^ ^ 

where An — w — t, A12 = A\ x = 0, A22 = uj — e, and B\\ 
is the Green function we are looking for. By working out 
the inverse matrix condition, we find in particular that 

AuBu + A 12 B 21 = 1 (F6) 
B21 = -A£A 21 B u (F7) 



and therefore 

(A u - A 12 A- 22 1 A 21 ) Bn = 1 . (F8) 
The Green function is thus 

G" 1 =Lj-t-9— — 0* 

UJ — £ 

= u-t- (F9) 
where we defined the so-called hybridization function 

r^H = £ °f°f (fio) 

Note that nowhere but in the last expression have we 
supposed that the matrix £ is diagonal. That condition 
simply serves to minimize computation time. 
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